Prepare libraies

suppressMessages(
suppressWarnings(
  c(library(data.table),
    library(tidyverse),
    library(rtracklayer),
    library(DESeq2),
    library(Biostrings),
    library(genomation),
    library(RColorBrewer),
    library(CLIPanalyze))
  )
)
##   [1] "data.table"           "stats"                "graphics"            
##   [4] "grDevices"            "utils"                "datasets"            
##   [7] "methods"              "base"                 "forcats"             
##  [10] "stringr"              "dplyr"                "purrr"               
##  [13] "readr"                "tidyr"                "tibble"              
##  [16] "ggplot2"              "tidyverse"            "data.table"          
##  [19] "stats"                "graphics"             "grDevices"           
##  [22] "utils"                "datasets"             "methods"             
##  [25] "base"                 "rtracklayer"          "GenomicRanges"       
##  [28] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
##  [31] "BiocGenerics"         "parallel"             "stats4"              
##  [34] "forcats"              "stringr"              "dplyr"               
##  [37] "purrr"                "readr"                "tidyr"               
##  [40] "tibble"               "ggplot2"              "tidyverse"           
##  [43] "data.table"           "stats"                "graphics"            
##  [46] "grDevices"            "utils"                "datasets"            
##  [49] "methods"              "base"                 "DESeq2"              
##  [52] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
##  [55] "Biobase"              "rtracklayer"          "GenomicRanges"       
##  [58] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
##  [61] "BiocGenerics"         "parallel"             "stats4"              
##  [64] "forcats"              "stringr"              "dplyr"               
##  [67] "purrr"                "readr"                "tidyr"               
##  [70] "tibble"               "ggplot2"              "tidyverse"           
##  [73] "data.table"           "stats"                "graphics"            
##  [76] "grDevices"            "utils"                "datasets"            
##  [79] "methods"              "base"                 "Biostrings"          
##  [82] "XVector"              "DESeq2"               "SummarizedExperiment"
##  [85] "DelayedArray"         "matrixStats"          "Biobase"             
##  [88] "rtracklayer"          "GenomicRanges"        "GenomeInfoDb"        
##  [91] "IRanges"              "S4Vectors"            "BiocGenerics"        
##  [94] "parallel"             "stats4"               "forcats"             
##  [97] "stringr"              "dplyr"                "purrr"               
## [100] "readr"                "tidyr"                "tibble"              
## [103] "ggplot2"              "tidyverse"            "data.table"          
## [106] "stats"                "graphics"             "grDevices"           
## [109] "utils"                "datasets"             "methods"             
## [112] "base"                 "genomation"           "grid"                
## [115] "Biostrings"           "XVector"              "DESeq2"              
## [118] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [121] "Biobase"              "rtracklayer"          "GenomicRanges"       
## [124] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
## [127] "BiocGenerics"         "parallel"             "stats4"              
## [130] "forcats"              "stringr"              "dplyr"               
## [133] "purrr"                "readr"                "tidyr"               
## [136] "tibble"               "ggplot2"              "tidyverse"           
## [139] "data.table"           "stats"                "graphics"            
## [142] "grDevices"            "utils"                "datasets"            
## [145] "methods"              "base"                 "RColorBrewer"        
## [148] "genomation"           "grid"                 "Biostrings"          
## [151] "XVector"              "DESeq2"               "SummarizedExperiment"
## [154] "DelayedArray"         "matrixStats"          "Biobase"             
## [157] "rtracklayer"          "GenomicRanges"        "GenomeInfoDb"        
## [160] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [163] "parallel"             "stats4"               "forcats"             
## [166] "stringr"              "dplyr"                "purrr"               
## [169] "readr"                "tidyr"                "tibble"              
## [172] "ggplot2"              "tidyverse"            "data.table"          
## [175] "stats"                "graphics"             "grDevices"           
## [178] "utils"                "datasets"             "methods"             
## [181] "base"                 "CLIPanalyze"          "GenomicAlignments"   
## [184] "Rsamtools"            "GenomicFeatures"      "AnnotationDbi"       
## [187] "plyr"                 "RColorBrewer"         "genomation"          
## [190] "grid"                 "Biostrings"           "XVector"             
## [193] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
## [196] "matrixStats"          "Biobase"              "rtracklayer"         
## [199] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [202] "S4Vectors"            "BiocGenerics"         "parallel"            
## [205] "stats4"               "forcats"              "stringr"             
## [208] "dplyr"                "purrr"                "readr"               
## [211] "tidyr"                "tibble"               "ggplot2"             
## [214] "tidyverse"            "data.table"           "stats"               
## [217] "graphics"             "grDevices"            "utils"               
## [220] "datasets"             "methods"              "base"

Load miR-peaks dataset

mirs.peaks <- readRDS("Datafiles/miRNA-peaks-list-12232019.rds")

Analysis

Prepare bw file for visualization

# combine bigwig files of the same genotype
$ bigWigMerge Exp-HF1_LIB044916_GEN00173026_R1_barcode.uniq.bw Exp-HF2_LIB044916_GEN00173029_R1_barcode.uniq.bw Exp-HF3_LIB044916_GEN00173032_R1_barcode.uniq.bw HF_merged.bedGraph

$ bigWigMerge Exp-HFK1_LIB044916_GEN00173035_R1_barcode.uniq.bw Exp-HFK2_LIB044916_GEN00173038_R1_barcode.uniq.bw Exp-HFK3_LIB044916_GEN00173041_R1_barcode.uniq.bw HFK_merged.bedGraph

convert chromosome name from Ensembl to UCSC format.

chromosomename <- read.table("../Joint_analysis/Ensembl.mm10.chrom.sizes.txt", sep = "\t")[,1][1:22]
hf_bed <- read.table("../Joint_analysis/HF_merged.bedGraph", sep = "\t")
hf_bed <- hf_bed %>% filter(V1 %in% chromosomename)
hf_bed$V1 <- paste0("chr", hf_bed$V1)
hf_bed$V1[hf_bed$V1 == "chrMT"] <- "chrM"
write_delim(hf_bed, file = "../Joint_analysis/HF_merged_edited.bedGraph", delim = "\t", col_names = F)

hfk_bed <- read.table("../Joint_analysis/HFK_merged.bedGraph", sep = "\t")
hfk_bed <- hfk_bed %>% filter(V1 %in% chromosomename)
hfk_bed$V1 <- paste0("chr", hfk_bed$V1)
hfk_bed$V1[hfk_bed$V1 == "chrMT"] <- "chrM"
write_delim(hfk_bed, file = "../Joint_analysis/HFK_merged_edited.bedGraph", delim = "\t", col_names = F)
# bedGraph to bigWig conversion
$ awk 'NR!=1' HF_merged_edited.bedGraph > input.HF_merged.bedGraph
$ sort -k1,1 -k2,2n input.HF_merged.bedGraph > sorted.input.HF_merged.bedGraph
$ bedGraphToBigWig sorted.input.HF_merged.bedGraph mm10.chrom.sizes HF_merged.bw

$ awk 'NR!=1' HFK_merged_edited.bedGraph > input.HFK_merged.bedGraph
$ sort -k1,1 -k2,2n input.HFK_merged.bedGraph > sorted.input.HFK_merged.bedGraph
$ bedGraphToBigWig sorted.input.HFK_merged.bedGraph mm10.chrom.sizes HFK_merged.bw
peaks <- mirs.peaks
peaks <- lapply(peaks, FUN = function(x) {x[order(x$count, decreasing = T)]})
peaks <- lapply(peaks, FUN = function(x) {x[x$padj<=0.01]})
mybw.dir <- "../Joint_analysis"
mybw.files <- list.files(mybw.dir, pattern = "bw$", full.names = T)[7:8]
print(mybw.files)
## [1] "../Joint_analysis/HF_merged.bw"  "../Joint_analysis/HFK_merged.bw"
cols <- brewer.pal(name = "Set2", n = 8)
reds <- brewer.pal(name = "Reds", n = 9)
blues <- brewer.pal(name = "Blues", n = 9)
mycolors <-c(blues[6], reds[6]) #HF, HFK
light.colors <- alpha(c(blues[3], reds[2]), 0.4)

Histogram

histogram plot function

#plot histograms
peaks_meta <- function(mypeaks = peaks, 
                       miRNA_family = "miR-451a", 
                       dispersion = "se",
                       dispersion.col = NULL,
                       coordinates = c(-400, 400), 
                       line.col = mycolors, 
                       winsorize = c(0,99),
                       title = ""){
  
  suppressWarnings(mypeaks <- lapply(mypeaks, FUN = function(x) {resize(x, width = sum(abs(coordinates)), fix="center")}))
  mypeaks <- GRangesList(mypeaks)
  mysml <- ScoreMatrixList(targets=mybw.files, window=mypeaks[[miRNA_family]], type = "bigWig", strand.aware = T)
  mysampleInfo <- data.frame(basename(mybw.files), c("HF", "HFK"))
  names(mysampleInfo) = c("sample", "genotype")
  names(mysml) = mysampleInfo$genotype[match(names(mysml), mysampleInfo$sample)]
  plotMeta(mysml, profile.names = c("HF", "HFK"), xcoords = coordinates, dispersion = dispersion, main =title, line.col = line.col, winsorize = winsorize, dispersion.col = dispersion.col)
}
len <- sapply(mirs.peaks, function(x) length(x))
mirs.peaks <- mirs.peaks[order(-len)]
mirna <- names(mirs.peaks)

for (i in 1:10) {
  peaks_meta(mypeaks = peaks, dispersion  = NULL, miRNA_family = mirna[i], title = mirna[i])
  peaks_meta(mypeaks = peaks, dispersion  = "se", miRNA_family = mirna[i], title = mirna[i],
           dispersion.col = light.colors)
}

Heatmap

plot heatmaps

peaks_heat <- function(mypeaks = peaks, 
                       miRNA_family = "miR-451", 
                       col = blues9, 
                       coordinates = c(-400, 400), 
                       order_rows = F, 
                       winsorize_parameters = c(1,98)){
  suppressWarnings(mypeaks <- lapply(mypeaks, FUN = function(x) {resize(x, width = sum(abs(coordinates)), fix="center")}))
  mypeaks <- GRangesList(mypeaks)
  mysml <- ScoreMatrixList(targets=mybw.files[c(1,2)], window=mypeaks[[miRNA_family]], type = "bigWig", strand.aware = T)
  mysampleInfo <- data.frame(basename(mybw.files), c("HF", "HFK"))
  names(mysampleInfo) = c("sample", "genotype")
  names(mysml) = paste(mysampleInfo$genotype[match(names(mysml), mysampleInfo$sample)], miRNA_family , sep = "_")
  mysml.scaled = scaleScoreMatrixList(mysml)
  #multiHeatMatrix(mysml.scaled, xcoords = coordinates, col = col)
  multiHeatMatrix(mysml, common.scale = T, xcoords = coordinates, winsorize = winsorize_parameters, col = col, order = order_rows)
}
colfunc <- colorRampPalette(c("white", "blue"))
mycols <- colfunc(128)
for (i in 1:10) {
  peaks_heat(mypeaks = peaks, miRNA_family = mirna[i], col = blues9, winsorize_parameters = c(0,99))
}

SessionInfo

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] CLIPanalyze_0.0.10          GenomicAlignments_1.24.0   
##  [3] Rsamtools_2.4.0             GenomicFeatures_1.40.1     
##  [5] AnnotationDbi_1.50.3        plyr_1.8.6                 
##  [7] RColorBrewer_1.1-2          genomation_1.20.0          
##  [9] Biostrings_2.56.0           XVector_0.28.0             
## [11] DESeq2_1.28.1               SummarizedExperiment_1.18.2
## [13] DelayedArray_0.14.1         matrixStats_0.58.0         
## [15] Biobase_2.48.0              rtracklayer_1.48.0         
## [17] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2        
## [19] IRanges_2.22.2              S4Vectors_0.26.1           
## [21] BiocGenerics_0.34.0         forcats_0.5.1              
## [23] stringr_1.4.0               dplyr_1.0.4                
## [25] purrr_0.3.4                 readr_1.4.0                
## [27] tidyr_1.1.2                 tibble_3.0.6               
## [29] ggplot2_3.3.3               tidyverse_1.3.0            
## [31] data.table_1.13.6          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0       ellipsis_0.3.1         fs_1.5.0              
##  [4] rstudioapi_0.13        farver_2.0.3           bit64_4.0.5           
##  [7] lubridate_1.7.9.2      xml2_1.3.2             splines_4.0.3         
## [10] cachem_1.0.1           impute_1.62.0          geneplotter_1.66.0    
## [13] knitr_1.31             jsonlite_1.7.2         seqPattern_1.20.0     
## [16] broom_0.7.4            gridBase_0.4-7         annotate_1.66.0       
## [19] dbplyr_2.0.0           compiler_4.0.3         httr_1.4.2            
## [22] biosignals_0.1.0       backports_1.2.1        assertthat_0.2.1      
## [25] Matrix_1.3-2           fastmap_1.1.0          cli_2.3.0             
## [28] htmltools_0.5.1.1      prettyunits_1.1.1      tools_4.0.3           
## [31] gtable_0.3.0           glue_1.4.2             GenomeInfoDbData_1.2.3
## [34] reshape2_1.4.4         rappdirs_0.3.3         Rcpp_1.0.6            
## [37] cellranger_1.1.0       vctrs_0.3.6            xfun_0.20             
## [40] rvest_0.3.6            lifecycle_0.2.0        XML_3.99-0.5          
## [43] zlibbioc_1.34.0        scales_1.1.1           BSgenome_1.56.0       
## [46] hms_1.0.0              curl_4.3               yaml_2.2.1            
## [49] memoise_2.0.0          biomaRt_2.44.4         stringi_1.5.3         
## [52] RSQLite_2.2.3          highr_0.8              genefilter_1.70.0     
## [55] plotrix_3.8-1          BiocParallel_1.22.0    rlang_0.4.10          
## [58] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [61] lattice_0.20-41        bit_4.0.4              tidyselect_1.1.0      
## [64] magrittr_2.0.1         R6_2.5.0               generics_0.1.0        
## [67] DBI_1.1.1              pillar_1.4.7           haven_2.3.1           
## [70] withr_2.4.1            survival_3.2-7         RCurl_1.98-1.2        
## [73] modelr_0.1.8           crayon_1.4.0           KernSmooth_2.23-18    
## [76] BiocFileCache_1.12.1   rmarkdown_2.6          progress_1.2.2        
## [79] locfit_1.5-9.4         readxl_1.3.1           blob_1.2.1            
## [82] reprex_1.0.0           digest_0.6.27          xtable_1.8-4          
## [85] openssl_1.4.3          munsell_0.5.0          askpass_1.1